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Abstract. - The square lattice with central forces between nearest neighbors is isostatic with a 
subextensive number of floppy modes. It can be made rigid by the random addition of next-nearest 
neighbor bonds. This constitutes a rigidity percolation transition which we study analytically by 
mapping it to a connectivity problem of two-colored random graphs. We derive an exact recurrence 
equation for the probability of having a rigid percolating cluster and solve it in the infinite volume 
limit. From this solution we obtain the rigidity threshold as a function of system size, and find 
that, in the thermodynamic limit, there is a mixed first-order-second-order rigidity percolation 
transition at the isostatic point. 
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Introduction. — Central force rigidity percolation de- 
scribes how a system of lattice sites can become rigid by 
randomly populating the bonds with springs or struts that 
can transmit central forces between sites. It features a 
mechanical critical point at which an infinite rigid cluster 
emerges and the system gains rigidity [l}|3]. Comparing 
to the analog scalar problem of connectivity percolation, 
in which an infinite connected cluster emerges and the 
system becomes conducting by populating bonds, rigid- 
ity percolation is a vector problem whose basic degrees of 
freedom are the d-dimensional position vectors of the lat- 
tice sites. Another important character of rigidity perco- 
lation is its long-range nature, that adding a bond at one 
place in the network could affect the rigidity of regions 
arbitrarily far away from that bond l2Jl3], which makes 
rigidity percolation a challenging problem for theoretical 
study. Numerical simulations on rigidity percolation re- 
vealed highly nontrivial physics. In two dimensions, nu- 
merical simulations using the "pebble game" algorithm on 
generic networks strongly suggest that the rigidity perco- 
lation is second order [2||4|. However the possibility of a 
weakly first order transition pip] is not completely ruled 
out due to finite size effects in the simulations. In three di- 
mensions, numerical simulations using the "pebble game" 
algorithm found rigidity percolation to be first order [71 . 

One special class of rigidity percolation occurs on pe- 
riodic isostatic lattices with the nearest-neighbor (NN) 
bonds already present from the start fsl . These lattices are 



at the onset of mechanical rigidity (isostaticity) because 
they have equal numbers of degrees of freedom and con- 
straints in the bulk, which for central forces corresponds to 
a coordination number z = 2d, where d is the dimension- 
ality of the system [o] . Examples are the two-dimensional 
square and kagome lattices and the three-dimensional cu- 
bic lattice. In a finite isostatic lattice, sites on the bound- 
ary have coordination number lower than 2d, which gives 
rise to a subextensive number of deformation modes in 
which none of the bonds change length. These so-called 
floppy modes do not cost any elastic energy, and are ex- 
tended across the lattice: the system cannot be macro- 
scopically rigid unless all floppy modes are somehow re- 
moved. To restore rigidity, one can randomly add some of 
the next-nearest-neighbor (NNN) bonds [4l[5[ [To}fT2] . This 
constitutes a special kind of rigidity percolation problem 
because all lattice sites are already connected and one only 
needs to remove a subextensive number of floppy modes. 
One therefore expects that a subextensive number of NNN 
bonds is enough to provide rigidity. 

In this Letter, we investigate rigidity percolation on the 
square lattice, which has coordination number z = 4 = 2d, 
and is isostatic [8 12 . We develop an analytical method 
to calculate the probability T{L,p) that a size L x L 
square lattice is made rigid by populating the NNN bonds 
with probability p, and thus we arrive at the rigidity 
threshold pn{L) that half of the realizations are rigid, i.e., 
F{L,p^{L)) = 1/2. Surprisingly, we find that the intu- 
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Fig. 1: (a) The variables Ai and Bi in which we express the zero energy deformations of the square lattice (before adding 
NNN bonds) are the relative displacements of entire neighboring rows and columns, respectively. The effect of a nonzero A2 
is illustrated, (b) Adding a NNN bond at location {i,j) corresponds to constraining the variables Ai and Bj to be equal, 
illustrated using B2 = A2 = B4 = x. (c) Graph representation of the network shown in (b), where each variable is represented 
by a node (light green and dark red for the Ai and Bi, respectively) and each NNN bond by a link between the corresponding 
nodes. 



itive argument used in earlier works plplnT], that having 
at least one NNN bond in each row and column should 
provide rigidity, presents a necessary but insufficient con- 
dition for rigidity of square lattices. In fact, one can add 
as few as L — 1 NNN bonds to an L x L square lattice to 
make the "one-per-row" condition satisfied, but according 
to the Maxwell counting ^ there are still L — 2 floppy 
modes. 

To obtain the rigidity threshold, we map the problem 
into a connectivity problem of two-colored (also known as 
bipartite) random graphs with L— 1 nodes of each color, for 
which the number of connected clusters corresponds to the 
number of floppy modes in the square lattice. We derive a 
recurrence relation for the probability that all nodes of the 
graph are connected into one cluster, which corresponds to 
the probability that the entire square lattice is rigid. We 
solve this equation asymptotically in the large L limit, and 
show that, (i) for nonzero p, F{L,p) ^ l — 2L{l — p)^, (ii) 
the rigidity threshold pn approaches zero when L — )■ 00 as 
Pr{L) ^ \nL/L + 0{1/L). We confirm these results using 
numerical simulations. 

Consider a regular square lattice with L sites per row or 
column. According to Maxwell ^, the number of floppy 
modes can be found by subtracting the number of con- 
straints from the number of degrees of freedom (usually, 
the global translations and rotations are subtracted as 
well, leaving only the nontrivial floppy modes). Thus the 
number of floppy modes in the L x L square lattice is 
2L^ ~ 2L{L — 1) — 3 — 2L — 3 which can be understood as 
follows. Each of the rows or columns of square plaquettes 
can be deformed into rhombi without changing the shape 
of the plaquettes in other rows or columns, as shown in 
fig. [T^. These deformations correspond to linearly inde- 
pendent floppy modes of the lattice. The total number of 
such rows and columns of plaquettes is 2L — 2. Note that 
any pure shear deformation can be obtained using linear 
combinations of these modes, which is why the shear mod- 
ulus is zero [^ Because the plaquette deformations concern 



relative displacements of neighboring rows or columns, the 
space they span does not contain the global translations 
(but global rotation is included). 

To remove these floppy modes and restore shear rigidity, 
we can add 2L — 3 NNN bonds to, for example, the left 
and the lower boundary plaquettes of the lattice. However, 
for the case of randomly populating each NNN bond with 
probability p, some of the NNN bonds could be redundant, 
and thus one may need more than 2L — 3 NNN bonds to 
restore the shear rigidity. Determining the rigidity of a dis- 
ordered network of springs always involves such counting 
arguments comparing the number of degrees of freedom 
to the number of independent constraints, and because 
of the disorder it is generally nontrivial to properly keep 
track of which constraints are truly independent, which 
makes rigidity percolation difficult theoretically [2|[3] . Be- 
low, we will derive simple rules for the removal of these 
floppy modes that make rigidity percolation on the square 
lattice tractable. 

The spring networks obtained through the addition of 
NNN bonds, while the NN bonds of the square lattice are 
already present from the start, are fundamentally different 
from the classical case of rigidity percolation on the tri- 
angular lattice, where no bonds are present when p = 0. 
Most strikingly, because all sites are already connected, 
the probability that a site belongs to the rigid cluster 
equals one as soon as p is large enough to have a span- 
ning rigid cluster. Thus, in the thermodynamic limit, the 
order parameter, defined as the probability that an arbi- 
trarily chosen site belongs to the percolating rigid cluster, 
jumps from to 1 at the transition, characterizing a first 
order transition. On the other hand, the system displays 
a diverging isostatic length scale I* ~ 1/p and a vanish- 
ing shear modulus fi ^ p^ [12] , which are characters of a 
second order transition. The discontinuous change in the 
order parameter and the continuous change in the shear 
modulus can be heuristically understood by realizing that 



^Due to anisotropy, strictly speaking, the square lattice has three 



elastic moduli, Cii, C12 and C44. With NN bonds only we have 
C12 = C'44 = 0, and in this paper we refer to C44 as shear modulus 
for convenience 1131. 
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the addition of a small number of bonds restores the rigid- 
ity of the whole system, thus these bonds carry all stress 
applied to the system, resulting in a vanishingly small 
shear modulus pi. This feature of a mixed first-ordcr- 
second-order phase transition resembles the jamming tran- 
sition in packings of frictionless spheres, which features a 
jump in the average coordination number from zero to the 
isostatic value, as well as power law scalings of the elastic 
moduli and the isostatic length scale 14,15 . The connec- 



tion between the square lattice and the jamming transition 
has been discussed in Refs. f8llT2l. 



The probability of rigid configurations. A 

square lattice of (m -I- 1) x (n -I- 1) sites with randomly 
added next nearest neighbor bonds is rigid if the only zero- 
energy modes are the global translations and rotations. To 
calculate the probability that this is the case, we consider 
the space of floppy modes of the original square lattice, 
which is a subspace of the space of all modes of the lat- 
tice. In particular, in linear elasticity the space of the 
floppy modes is the null space of the dynamical matrix of 
the lattice. To exclude the trivial global translational de- 
grees of freedom, we choose to project this space of floppy 
modes into the basis of deformations of the square plaque- 
ttes in each row or column into rhombi, as discussed earlier 
and shown in fig. [I^, or more precisely, the relative hori- 
zontal (vertical) displacements between neighboring rows 
(columns). The dimension of this space is m -I- n with 
the global rotation included. This leads to a description 
in terms of relative row-displacements Ai, A2, . . . , Am and 
relative column-displacements _Bi, i?2, • . • , Bn- 

As shown in fig. Ilk, exciting the floppy mode labelled by 
A2 turns all squares on the corresponding lattice row into 
rhombi. To make the system rigid, all floppy modes have 
to be constrained, so that one needs to have at least one 
NNN bond on each row and on each column. However, as 
will soon become clear, this necessary condition for rigidity 
is not sufficient. 

Adding a NNN bond to the square at (i, j) constrains 
it to remain square in any floppy deformation: it is only 
allowed to rotate. This amounts to setting Ai = Bj (with 
the proper choice of signs). Each NNN bond provides 
such an equality, as illustrated in fig. [TJd, in which two 
NNN bonds set B2 = ^2 = -84 = x. The rigidity of the 
whole lattice would correspond to having Ai = Bj — x for 
alH = 1, . . . , TO and j — 1, . . . ,n, with x corresponding to 
the amount of global rotation in this case. This suggests 
a graph representation of realizations of the rigidity per- 
colation procedure. The nodes represent the variables Ai 
and Bj, and edges between two nodes represent equalities 
of the corresponding variables that arise from the NNN 
bonds (fig. [IJ:). Hence, we obtain a description of rigidity 
percolation in terms of a random graph with two types of 
nodes (green A nodes and red B nodes), with edges be- 
tween an A-node and a _B-node present with probability 
p. This mapping has been employed before in the context 
of studying rigidity in square and cubic structures with 



added NNN bonds, albeit without the probabilistic aspect 
of rigidity percolation 
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Within this description, a rigid lattice is represented by 
a connected graph, i.e., a graph in which there is a path 
along edges between any given pair of nodes, correspond- 
ing to the case of all Ai = Bj with global rotation left 
as the only degree of freedom. Note that we use cross- 
bars to represent the effect of the NNN bonds, because 
the second NNN bond on a plaquette is redundant, when 
only the question "rigid or not" is considered. Near the 
percolation transition in large systems where p^ — >■ 0, it 
is equivalent to use crossbars with probability p or to use 
separate NNN bonds with probability p/2, so the resulting 
scalings are the same. 

We now derive the scaling of the probability T{m, n,p) 
that the system is rigid, as rn, n — >■ 00. This can be ob- 
tained as a corollary from the results of Palasti [Tt], but 
we choose to present our own, more accessible derivation. 



For the usual Erdos-Renyi model 18 ■ 20 for generating 



random graphs in which there are n equivalent nodes and 
each of the n{n — l)/2 edges connecting any two nodes is 
present with probability p, the probability Ti{n,p) that 
the graph is connected can be calculated recursively from 
an expression due to Gilbert 
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1- J"i(n,p) = ^ 



fc=i 



J'i{k,p)q 



k{n—k) 



(1) 



where q = 1 — p. Here, the probability that the graph is 
not connected is written as the sum over all possible sizes 
k of the cluster that an arbitrarily chosen node (labelled as 
node 1) is part of. The binomial coefficient represents the 
number of ways to choose the other nodes in the cluster, 
and the power of q denotes the probability that none of 
those k nodes are connected to any of the other n—k nodes. 
Gilbert showed that for large n, Fi{n,p) 
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l-nq" 

so the connectivity threshold for the usual Erdos-Renyi 
model is pc{L) ^ \nn/n. 

We adapt eq. (fTl) for our two-colored graph in which 
edges between differently colored nodes are present with 
probability p. Again, we consider the sum over all the 
possibilities for the size of the cluster having k green nodes 
and / red nodes that an arbitrary green node (to which we 
assign the variable ^1) is part of: 
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A recurrence relation for T{m,n,p) can be obtained by 
moving the term for k = m,l = n from the sum to the left 
hand side of the equation. In the Appendix we show that 
as TO = n — > cx) the probability approaches 



J-{n, n,p) sa 1 — 2nq^' 



00 



(3) 



For any finite p the probability approaches unity, and the 
value of p needed to make half the realizations rigid ap- 
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Fig. 2: (a) Numerically obtained fraction of rigid configurations 
J-{n, n, p) as a function of bond-presence probability p for var- 
ious system sizes. The solid lines denote the asymptotic form 
T = 1 — 2n(l — p)" [eq. ([3|]. For completeness we included a 
few points (open triangles) obtained by numerically evaluating 
the recurrence relation S for n — 80. (b) The same data with 
p rescaled by ln(2n)/n, according to eq. Q. 



proaches zero as (see Appendix) 
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+ 0{l/n) as n — > oo . 
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Numerics. — For various values of n and p, we gen- 
erated 1000 trials of the rigidity percolation procedure 
and analyzed the corresponding graphs using an adapted 
version of the Hoshen-Kopelman cluster labelling algo- 
rithm 
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We plot the fraction J^ of the graphs that 
are connected, corresponding to the fraction of rigid net- 
works, in fig. [2^. The solid lines in this figure indicate 
the asymptotic form derived in eq. ([3]) . Figure Pp shows 
the same data, with the p-axis rescaled by ln(2n)/n [the 
2 comes from the minimum prefactor of the 1/n term in 
eq. Q — see Appendix], clearly showing the scaling of 
Pnin). 

Discussion. — The graph picture for rigidity percola- 
tion on the square lattice provides a transparent way to 
keep track of the effect of added constraints. The intu- 




Fig. 3: Having a crossbar on each row and each column does not 
guarantee rigidity, (a-e) The five floppy modes of an example 
system, (f ) The graph representation of this example system. 
The five connected clusters correspond to the five floppy modes 
in panels a-e of this figure. 



at least one NNN bond in each row and column should 
provide rigidity, corresponds to not having any isolated 
nodes with no edges. Clearly this is a weaker condition 
than demanding the graph to be connected. This is illus- 
trated in fig. |3J in which panels (a-e) show the five floppy 
modes of a lattice in which the condition of having at least 
one NNN bond in each row and column is satisfied. The 
corresponding graph is shown in fig. Isf. While counting 
the number of floppy modes by just looking at a picture of 
a lattice is far from trivial, the graph representation gives 
a much simpler view in which it is easy to see that there 
are 5 distinct connected components, so that the lattice 
has 5 floppy modes. 

Thus, using the weaker "one-per-row" condition leads to 
a structural underestimation of the rigidity threshold pR. 
However, this condition does provide the correct scaling 
with system size, pR 



hiL/L 11 . This implies that 



the fraction of realizations that satisfy the "one-per-row" 
condition but are not fully connected vanishes as L -^ 
oo, which is consistent with the known result that just 
below the transition, the graph contains just one giant 
component and some isolated nodes 22 . 

As an aside, we note that the appearance of this gi- 
ant component is in itself an important transition in the 
theory of random graphs 
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itive argument used in earlier works (4l5 11 , that having 



However, for the system 
studied in this paper its relevance is limited. Having a 
giant component is not enough for rigidity to percolate, 
because any isolated node that is not connected to the 
giant component corresponds to a row or column in the 
lattice in which there are no crossbars, so that the two 
parts of the system on either side of that row or column 
can freely shear past each other. Thus the relevant tran- 
sition is indeed that to connected graphs, as described in 
this letter. 
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The formulation of our approach in terms of the null 
space of the dynamical matrix implies that the floppy 
modes we talk about in principle only need to be floppy up 
to linear order, i.e. for infinitesimal deformations. How- 
ever, the row and column displacements in our discussion 
can be straightforwardly extended to finite displacements, 
and thus the mapping to the graph representation is not 
limited to linear elasticity. This can be easily seen in fig. |3J 
in which the modes have a finite amplitude but are clearly 
zero energy modes. 

In our discussion the square lattice has been assumed to 
be a perfect periodic lattice. However, "generic lattices" , 
in which only the connectivity (topology) is prescribed but 
bond lengths are not all equal, are often used in the dis- 
cussion of rigidity theory [3||23] . The percolation problem 
is still well-defined: whether or not a given structure (a 
set of points and rigid bars connecting them) is rigid does 
not depend on the precise positions of the points, but only 
on the connectivity and on the fact that the points do not 
have any positional order [^ [23] . This, again, seems to 
suggest some kind of graph theoretical approach could be 
helpful. However, the mapping to the graph connectivity 
problem that we used in this paper relies specifically on 
the positional order of the sites and hence does not work 
for the case of generic lattices. The reason is that some of 
the clusters of the resulting graph represent internally in- 
consistent constraints in this language. These clusters do 
not correspond to fioppy modes, so that the lattice could 
become rigid before the graph is fully connected. The 
threshold probability pB.{n) for periodic lattices therefore 
serves as an upper bound for the pr{ti) for generic lattices, 
but there is no obvious lower bound, because the "one-per- 
row" condition that provided the lower bound on pR in the 
perfect square lattice (see appendix) is no longer a neces- 
sary condition for rigidity. Clearly, a more sophisticated 
approach is needed to describe generic rigidity percolation. 
The opportunities that arise from mapping this problem 
to so-called hypergraphs are currently being explored j24 j. 

In conclusion, we have derived an exact recurrence equa- 
tion for the probability jF(L,L,p) that a square lattice 
with given size L and NNN bonds occupation probabil- 
ity p is rigid, and obtained the asymptotic solution to 
this equation in the limit of i — > oo. Our results unam- 
biguously prove that the rigidity percolation in the square 
lattice occurs at p = and is a mixed first-order-second- 
order transition: In the thermodynamic limit, on the one 
hand all lattice sites immediately become part of a perco- 
lating rigid cluster as soon as p is nonzero; on the other 
hand, the shear modulus increases continuously from zero 
at the transition ^. Our result indicates that the sys- 
tem shows no fractal spatial dimension in the percolating 
cluster. Therefore a mean-field approach, such as the co- 
herent potential approximation used in Ref. 12 , may be 



sufficient for the square lattice near its rigidity threshold. 



Our results also verify the existence of a diverging length 
scale I* ^ p~^, ignoring the slowly varying In factor. For 
any given NNN bonds occupation probability p, lattices 
of sizes bigger than I* are very likely to be rigid and lat- 
tices smaller than I* are very likely to be floppy. This 
observation in square lattices is consistent with the cut- 



ting argument on the isostatic length scale by Wyart 25 



Counting independent constraints is key in various sys- 
tems that show a rigidity transition, from jammed sphere 
packings 
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to rigidity percolation [l||3]. Mapping the 
rigidity problem into the connectivity problem of random 
graphs allows to draw inspiration from the vast body of 
work on (random) graphs to gain insight in the rigidity 
problem. This approach has led us to an exact expression 
for the probability that a square lattice is made rigid by 
randomly adding next nearest neighbor bonds. We spec- 
ulate that this idea can be applied to a wider range of 
models for random media, and advance our understand- 
ing of disordered materials. 
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Appendix: The infinite system size limit. To 

obtain the thermodynamic limit m,n —^ cx), we follow 
Gilbert's strategy of deriving a lower and upper bound 
to J-{n,n,p) and showing that these are in close agree- 
ment [18] . The upper bound on J-{n, n,p) is given by the 
lower bound on 1 —J-(n, n,p) that is set by the probability 
that at least one of the 2n nodes is not connected to any 
other node. Denoting by Ei the event that node i is not 
connected to any other node, which has probability g", we 
use a Bonferroni inequality [26' to obtain 

1-T{n,n,p) > P (U^'J > E^(^')-E^(^''^^j) ■ 

\ i / i i<j 

(5) 
The first term on the right hand side equals 2nq". The 
second term is of order (nq")^ and can therefore be ignored 
as n -^- oo, so that the upper bound becomes 



■F{n, n,p) < 1 ~ 2nq" as n — > oo 



(6) 



■^Genericness of a set of site positions is usually defined as there 
not being any nontrivial algebraic equation (with rational coeffi- 
cients) relating the coordinates of the sites. 



The lower bound is obtained directly from the recur- 
rence relation (l2|. We have, using J^(fc, l,p) < I and writ- 
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where in the first line the prime indicates the term (fc, /) — 
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higher. Hence we find convergence of this lower bound 
with the upper bound in eq. ([6|, and conclude that 
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This asymptotic form holds for finite p and tells us how 
T aproaches unity as bigger systems are considered. Note 
that the scaling of pjx., defined as the value of p where 
J^(n,n,p) = 1/2, does not immediately follow form this 
because p^ vanishes as n ^>- oo. However, we can obtain 
bounds on p^ by considering where the bounds on J^ cross 
1/2, as long as we evaluate them keeping all orders of ng". 
Equating the right hand side of eq. (l5| to 1/2 and solving 
for very large n we obtain q" — l/(2n), which gives a lower 
bound on pj^ of 



PR > 



ln2n 



(X) 



(9) 



The upper bound on p^ follows from equating the right 
hand side of eq. ^\ to 1/2. Numerically, we find again 
a solution of the form g" — l/(/3n), but this time with 
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From the result of the more elaborate proof by Palasti [17] 
one can derive that 



PR 



ln(2n/ln2) ln2.89n 



which falls nicely within the bounds we derived by simpler 
means. In any case the leading order behavior is given by 
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